Ionic and electronic structure of sodium clusters up to N=59 
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We determined the ionic and electronic structure of sodium 
clusters with even electron numbers and 2 to 59 atoms in ax- 
ially averaged and three-dimensional density functional cal- 
culations. A local, phenomenological pseudopotential that 
reproduces important bulk and atomic properties and facil- 
itates structure calculations has been developed. Photoab- 
sorption spectra have been calculated for Na2, Nag, and Na^ 
to Na^g. The consistent inclusion of ionic structure consid- 
erably improves agreement with experiment. An icosahedral 
growth pattern is observed for Na]''g to Na^g. This finding is 
supported by photoabsorption data. 

PACS: 36.40.Vz,31.15.Ew,71.15.H,36.40.Wa 



I. INTRODUCTION 

Since the pioneering experiments of Knight et 
and their interpretation in terms of the jellium model, 
sodium clusters have attracted great attention, both ex- 
perimentally and theoretically. This is due to the fact 
that sodium is the "simple metal" par excellence, and 
thus is best suited for the study of fundamental effects. 
From photoabsorption experiments it is known that 
small sodium clusters have overall shapes that strongly 
vary wftth size and are determined by electronic shell 
effects.Q eI On the other hand, cold clusters with several 
thousand atoms build icosahedraja i.e., they show ionic 
shell effects, whereas the bulk material crystallizes in a 
bcc lattice. Recent experiments indicate that both ionic 
and electronic degrees of freedom play a role in deter- 
mining the structural and jtliermal properties of clusters 
with several tens of atoms.t2l 

On the theoretical side, there have been different ap- 
proaches to obtain an understanding of the delicate inter- 
play between ions and electrons that is the source of this 
variety. On one hand, relatively transparent models like 
the shell and jellium models m.^sei/eml levels of sophistica- 
tion have been widely used.al3'ElrE£l On the other hand, 
quantum chemistry and density functional theory offer 
methods to study clusters on the highest level of sophis- 
tication presently possible. But ab initio calculations in 
the strictest sense, i.e., taking all electrons into account, 
have only been performed for a few selected cases for the 
smallest clusters, due to t^e enormous computational ef- 
fort that they require.Oilj Taking only the valence elec- 
trons into account reduces the complexity considerably. 



eaxching low-energy configurations 
~Ej still grows rapidly with the sys- 



but the expense foii^ 
in three dimensionsEB 
tem size. The largest ah initio studies of sodium clusters 
to our knowledge are the recent finite-temperature inves- 
tigations presented in Refs. |25| , ^ . 

Several models have been developed to bridge the gap 
between the ah initio calculations and the shell and jel- 
lium models. The "Spherically Averaged Pseudopoten- 
tial Scheme" (SAPS) describes the ions by pseudopoten- 
tials, in most cases local ones, and thepralence electrons 
are restricted to spherical symmetry.EZi Models based 
on a volui ne-iav eraged or perturbational treatment of 
ionic effectsE§Ej improve on the treatment of the elec- 
trons. Yet fuctiier approaches are the Hiickel and re- 
lated models InrH molecular dynamics based on empiri- 



cal potentials £3 and recently, also the extended Thomas- 
Fermi model combined with a localj^eudopotential has 
been used to study sodium clusters .EjO 

From the above examples, two points become clear 
that considerably facilitate systematical studies of clus- 
ters with more than twenty atoms. First, the fact that 
many calculations make use of phenomenological local 
pseudopotentials shows the importance of such poten- 
tials. This is jCspccially true for sodium, where previous 
investigationg|i30'Lj have shown that a local electron-ion 
interaction can be a good approximation. But especially 
in cluster physics, where one bridges the region between 
the atom and the bulk, it is important that a pseudopo- 
tential give reliable results for all sizes despite its locality. 
Second, one needs models that make it computationally 
manageable to calculate structures and optical properties 
of clusters systematically for a wide range of considerable 
sizes, but which on the other hand leave the underlying 
physics intact and are detailed enough so that relevant 
information can be drawn from them. Besides extending 
the computational range, such models will serve the even 
more important purpose to distill the dominant physical 
effects from the wealth of details that fully ab initio cal- 
culations supply. The results of the model calculations 
of course must be verified in calculations of higher accu- 
racy and by comparison with experiment. These tasks 
are addressed in the present article. 

In Section || we develop a new phenomenological pseu- 
dopotential which meets the just mentioned require- 
ments. Detailed comparisons with ab initio calculations 
for the smallest clusters in section [II verify the validity 



of the pseudopotential and an axially averaged density 
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functional mode]0 for structure calculations. In section 
IV , a systematic survey over structures and photoabsorp- 
tion spectra of Na clusters up to Na^g is given. The re- 
sults of the survey are summarized and discussed in the 
concluding section 0. 



II. A LOCAL PSEUDOPOTENTIAL FOR 
STRUCTURE CALCULATIONS 

Rigorous, pseudopotentials in the sense of PhillipSi-and. 
KleinmarH and modern ab initio pseudopotentiala3'c2l 
are always nonlocal. However, it has been noted early 
in the development of pseudopotential theory that by 
relaxing the Phillips-Kleinman condition, one can open 
up a new class of pseudopotentials .El These are also 
termed "model potentials" because they are constructed 
by choosing some analytical functions as models for the 
partial-wave potentials and adjusting their parameters 
such that some experimentally known quantities, e.g., 
atomic energy levels, are matched. If several partial- 
wave potentials can be chosen to be the same, one can 
construct a local model potential. In fact, such local phe- 
nomenological potentials have been and are being used 
successfully in many branches of physics, see, e.g., Refs. 
p9| , ^ , ^ , ^2|j4^ . However, the local approach can only be 
expected to be good for the so-called "nearly free elec- 
tron" metals whose atomic structure consists of only s or 
p electrons outside of an ionic core with a noble gas con- 
figuration. (Lithium, with its lack of p electrons in the 
core, is the prominent exception.) In these metals, the 
two contributions making up the pseudopotential have a 
tendency to cancel each other, and the effective poten- 
tial a valence electron experiences is further diminished 
by screening effects. Therefore, model potentials have 
been proposed in the past which tried to exploit this 
simple electron icj s.t|rijrt|Ure either by fitting pmperties of 
the bulk solidj30'OE3 or of the single atcmruo or by ad 
hoc fulfilling desired numerical properties E3'C3 However, 
some of these potentials can lead to wrong predictions 
when they are used in physical surroymdings different 
from the one in which they were set up,Ej or when prop- 
erties other than the adjusted ones are looked at. With 
an emphasis on solid state properties, the last point has 
been discussed in detail previously, see, e.g., Ref. 35 for 
an overview. Our aim here is to develop a local model 
potential that gives a maximum degree of transferability 
in the sense that potentials constructed according to our 
scheme should reproduce the important physical proper- 
ties of a system, irrespectively of its number of atoms or 
the way in which these are arranged. 



A. The Ansatz 

The construction of a model potential consists of two 
steps. The first is to choose the model function. It should 



allow for a physical interpretation of the final potential, 
and at the same time should have analytical and numer- 
ical properties that allow for an easy application. The 
parameterization 



Vps (r) 



where 



r 



C2erf 



V2a2 



2 r 

erf(a;) = / dy exp 



(2.1) 



(2.2) 



certainly meets the second requirement, since the errpi; 
function can very efficiently be handled numericalljicZI 
and yields a smooth representation on a grid. That this 
parameterization also has a transparent physical inter- 
pretation will be demonstrated at the end of this sec- 
tion. The second step is the choice of the four parameters 
CTi , (72 , ci , C2 . One necessary condition is that the correct 
Coulomb limit. 



hm Vps (r) = 

r — >oo 

be obtained, which requires 

C1+C2 ^ 1. 



(2.3) 



(2.4) 



Thus we are left with three free parameters. Since our 
aim is to construct a pseudopotential that gives reliable 
properties for all clusters, i.e., spanning the region from 
the atom to the bulk, the most important properties of 
both atom and bulk solid must be reproduced. There- 
fore, we have chosen to fit the Wigner-Seitz radius and 
the compressibility B of the crystalline metal, together 
with the energy of the atomic 3s level e^- These quanti- 
ties characterize sodium and influence both structure and 
electronic excitations, thus being of great importance for 
reliable results. A test for whether we really have cap- 
tured the relevant physics will be to check if non-fitted 
quantities (bond- lengths, atomic spectra, bulk binding, 
dipole resonances) are also reproduced correctly. 



B. Determination of basic properties 

We calculate the bulk properties in second-order per- 
turbation theoryHj The unperturbed system is the nonin- 
ter acting homogeneous gas of valence electrons, the per- 
turbation is given by the potentials arising from the crys- 
tal lattice and from the interaction of the electrons with 
each other. Each crystal ion is described by a pseudopo- 
tential Vps centered on a lattice site. Up to second order, 
the binding energy per valence electron eb is given by 



eb(?'s) = ekinlj's) + e^cirs) + epsi(rs) -I- e^irs) + en 



sirs). 



(2.5) 
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Here Ts = [3/(47rn)]^ is defined in terms of the average 
valence electron density n, 



ekin(?'s) 



10m 



1 



(2.6) 



is the noninteracting kinetic energy, and exc(?'s) is the 
exchange and correlation energy for which we haiie em- 
ployed the LDA functional of Perdew and Wang.H The 
quantity 



(rs) 



(2.7) 



is the first-order contribution of the pseudopotential to 
the binding energy, eh(?'s) is the Hartree energy of the 
valence electron density, and en is the electrostatic energy 
of point ions. Finally, ebs(?'s) is the band structure energy 
discussed below. Integrations are taken over the (infinite) 
volume of the crystal. These contributions are rewritten 
as 



epsi(?'s) + eh{rs) + eii(rs) 
3 f fVh ^e2\ 



(2.8) 



to obtain the volume-averaged repulsive part of the pseu- 
dopotential that defines its strength 



Ks + — UV, 



(2.9) 



and the total electrostatic energy of point ions in a 
compensating uniform negative background, called the 
Madelung energy e^. Separating thej-Coulomb force 
into long and short-range components,c2l one can cal- 
culate em, and for the bcc lattice one obtains., em = 



-0.895929^3 e^/rs. Thus the only characteristidSU of the 
pseudopotential that enters the binding energy in first- 
order perturbation theory is its strength Sps- For fixing 
the radial dependence of the potential it is therefore es- 
sential to take into account the second-order band struc- 
ture energy 



ehs{rs) 



1 



2 Z^Airr^ 



J2\Vps (q)5(q) 



In ( |2iq ). 



q#0 



N 



x(q) 



(2.10) 



•^(q) = ^I]exp(-zqR^- 



(2.11) 



is the structure factor with the sum running over all ionic 
positions. For a lattice without basis one obtains S (q) = 
(^q.K) where K denotes reciprocal lattice vectors, i.e., the 



sum in (2.10) is running over reciprocal lattice vectors 
only. Furthermore, 



X(q) 



27r2fi2 



27? 



l + V 



1-7? 



2kF 

(2.12) 



is the Lindhard function with hp = {3TT'^n)3 , and 



e(q) = i- 



(i-e(q))x(q) 



(2.13) 



is the dielectric function including the local field correc- 
tion 



47re2 dn'^ 



(2.14) 



in the LDA. The pseudopotential enters via its Fourier 
transform, which for our potential is given by 



Vps (q) 



■ ^ Ci exp 



(2.15) 



From Eqs. (U) to ( ^.13| ) we calculate eb, the minimum 
of eb determines rg, and the bulk modulus is obtained 
from 

dP 1 /a^eb 2 9eb\ 

The reciprocal lattice vectors are generated numerically 
from the reciprocal basis, and we carefully checked that 
in all calculations the sum over reciprocal lattice vectors 
was numerically converged. 

The atomic 3s level is the lowest eigenvalue ea of the 
Schrodinger equation 



^ Ijl + l) 



2m 
If 



Vp. (r) 



27n 



ea w (r) = 



(2.17) 



for I = 0, where the pseudo wavefunction has, as 
usual, been factorized into radial and angular compo- 
nents, ip{r) = [u{r)/r] Yi^{^,ip). At this point, a sub- 
tlety should be considered. The atomic energy calcu- 
lated within density functional theory using the LDA 
will slightly differ from ea due to the well known fact 
that in the LDA, the self interactions in the Hartree and 
exchange energy do not cancel each othjsr exactly and 
leave a spurious self interaction energy.Eil One therefore 
might be tempted to "compensate" this self-interaction 
energy by adjusting the pseudopotential parameters such 
that the experimental value is matched when the self- 
interaction energy is included. This, however, would be 
dangerous for several reasons. First, it must be recalled 
that for an accurate description of bonding properties, 
the pseudopotential must lead to an accurate description 
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of the electronic density. If one adjusts the pseudopoten- 
tial parameters such that a spurious energy is compen- 
sated, then the valence electron density resulting from 
this pseudopotential might correspondingly show spuri- 
ous, unphysical deformations, leading Juo wrong bond- 
ing properties, as discussed previously.E3 Secondly, the 
self-interaction energy becomes smaller with increasing 
delocalization of the electrons. For an electron delocal- 
ized Qj£fir a volume 57, the self-interaction correction is of 
orderEj Since the valence electrons in metal clus- 

ters are delocalized to a high degree, and since the cluster 
volume changes noticeably in the size range from N=8 
to N=58 that we are interested in, a compensation of 
the atomic self interaction energy via the pseudopoten- 
tial would lead to problems for increasing cluster sizes. 
And third, building the self-interaction energy into the 
pseudopotential would be inconsistent with the proce- 
dure of fitting to the bulk, because the bulk calculation 
is less affected by the self-interaction error due to the de- 
localization of the b ulk e lectrons. Therefore, it is better 
to exactly solve Eq. ( ^.17 ). This can straightforwardly be 
done numerically by combining the Runge-Kutta method 
with adaptive stepsize control with a globally convergent 
Newton scheme. As a welcomed side effect, our pseu- 
dopotential might also be usable in self-interaction free, 
i.e., beyond-LDA density functional calculations. 



C. Results and comparison 

Th e three free parameters of our model potential 
(2.1) were now chpsen such that the experimental low- 



temperature valueo = 3.93ao be reproduced exactly, 
while at the same time ea and—B be as close as possL 
ble to the experimental valuescil = —0.378Ry andEj 
B — 0.073AIbar. These conditions lead to the parame- 
ters 



cri= 0.681 ao, ci = -2.292, 
cr2 = 1.163ao, C2 = 3.292. 



(2.18) 



In Table | we have listed the resulting Vg, B and Ca 
for our smooth-core pseudopotential and for other local 
pseudopotentials that have widely been used in cluster 
physics. The empty-core potential with both of the most 
frequently used choices for its cut-off radius Tc, and the 
pseudopotential of Ref. ^ that was constructed in the 
same spirit and adjusted in first-order perturbation the- 
ory only, lead to considerable deviations from the ex- 
perimental values f(p,|aJ.]|-quantitics. The local Heine- 
Abarenkov potentialc3c3S gives a reasonable rg and an 
Ca very close to the experimental one, but 10% error in 
B. Our pseudopotential by construction gives no or only 
very small differences from the experimental values for Vg , 
B and ea, showing that it is possible to obtain reason- 
able results for all these quantities with a local potential. 
However, a severe test will be whether also non- fitted 
quantities are reproduced correctly. To this end, we have 



Pseudopotential 


Ts/ao 


B/Mbar 


ea/eV 


Empty core, rc= 1.66" 


3.49 


0.119 


-5.52 


Empty core, rc= 1.76'' 


3.61 


0.109 


-5.32 


Ref. il 


3.84 


0.079 


-5.31 


Heine- Abarenkov'^ 


3.90 


0.080 


-5.12 


Present work 


3.93 


0.074 


-5.18 


Experiment 


3.93 


0.073 


-5.14 



TABLE I. Equilibrium density rs, bulk modulus B, and 
atomic 3s level ea, calculated for standard local pseudopoten- 
tials: (a) Ref. 0, (b) Ref. and (c) Ref. ^|^||. Bottom 
line: Experimental results. See text for details and references 
for experimental values. 



calculated the bulk binding energy Cb, the dimer bind- 
ing length ddimcr, the zero of the pseudopotential form 
factor, given by go = ^2[ln(c2) — ln(— ci)]/ (cr| — crj) in 
our case, and the first seven excited atomic energy levels. 
Table |l| shows that also these ten non-fitted quantities 
come out very close to the measure d v alu es, r evealing 
that the model potential defined by (2T), ( ^.18| ) indeed 
incorporates the relevant physical effects. 





Present work 


Experiment 




-6.19 


-6.25 


go/(2fcF) 


0.92 


0.87°/0.98'' 


ddimcr/fflO 


5.78 


5.82 


3s - 3p /eV 


2.19 


2.10 


3s - 4s /eV 


3.18 


3.19 


3s - 3d /eV 


3.68 


3.52 


3s - 4p /eV 


3.80 


3.75 


3s - 5s /eV 


4.13 


4.12 


3s - 4d /eV 


4.33 


4.29 


3s - 5p /eV 


4.39 


4.35 



TABLE II. Left column: Bulk binding energy eb from per- 
turbation theory; normalized zero of the form factor go/(2fcF); 
dimer binding length ddimor calculated with our smooth-core 
pseudopotential and CAPS; valence electron excitation ener- 
gies for our pseudopotential. Right column: Measured bulk 
binding energy;Ej semi-empirical values for go / {^kp) from op- 
tical properties (a) and from elastic ccpstants (b), see Ref. ^ 
for a discussion; dimer binding length^a and psergy level dif- 
ferences as obtained from spectroscopic lines £j 

Fig. |l] shows our pseudopotential in real and in Fourier 
space, leading us back to the ques tion of the motivation 
for the ansatz (2.1). From Eq. (2.4), one can interpret the 
coefficients ci and C2 as the charges associated with the 
attractive and repulsive terms, respectively, in the pseu- 
dopotential. This interpretation becomes more transpar- 
ent when one looks at the corresponding pseudocharge 



density. The pseudocharge rips is related to the pseu- 
dopotential via Poisson's equation 



'ps 



(r) 



(2.19) 



and for our pseudopotential is given by 
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FIG. 1. The smooth-core pseudopotential in real and 
Fourier space. 



ni exp 



2a 



n2 exp 



where 



for 



1,2. 



, (2.20) 



(2.21) 



Thus, our parameterization describes the ionic core by 
two overlapping Gaussian charge densities. The different 
heights and widths of the Gaussians result in a pseudo- 
density that is negative for small distances, correspond- 
ing to a repulsive core, and has a positive tail that com- 
pensates the core for larger distances. Equation ( 2.20| ) 
thus is a generalization of thej)arameterization that was 
used in our previous studies Jl3 where the pseudodensity 
has a two step-profile. Besides its tr ans parent physical in- 
terpretation, the pseudopotential (2.1) has excellent nu- 
merical properties: in real space the potential is smooth 
with no extreme peaks and can very efficiently be han- 
dled numerically via the pseudodensity, allowing to solve 
the Coulomb problem for electrons and ions together. At 
the same time, its rapidly converging Fourier transform 
makes the potential equally well suited for calculations 
in reciprocal space. 



III. COMPARISON WITH AB INITIO RESULTS 

In the present work we have employed the "Cylindri- 
cally Averaged Pseudopotential Scheme" (CAPS), which 



has been introduced earlieri3@ and which has been mod- 
ified and improved for the present purposes. Since CAPS 
treats the electronic degrees of freedom in the Kohn- 
Sham formalism without limiting them to spherical sym- 
metry, explicitly includes the ionic structure, and fur- 
thermore is numerically efficient, it is a very good tool 
for systematical studies of the interplay between ions and 
electrons. But before CAPS and the smooth-core pseu- 
dopotential are applied on a large scale, we further test 
their usefulness for cluster structure calculations. 

That the smooth core pseudopotential reproduces the 
experimentally known dimer bond length has already 
been shown in Table ||. For the next larger clusters, no 
direct experimental information on structures or bond 
lengths is available, and we therefore resort to compar- 
isons with other theoretical work. Fig. ^ shows the four 
smallest Na clusters with even electron numbers as they 
are obtained in CAPS with the smooth core pseudopo- 
tential. For these clusters, also three-dimensional geom- 
etry optimizations have been performed using ah initio 
pseudopotentiala apd Hartree-Fock with configuration in- 
teraction (CI)BB, or DFT with the LDAB^EI for the 
valence electrons, respectively. Also all-jelectron Hartree- 
Fock calculations have been reported.EZlL3 CAPS finds 
exactly the same structures as the three-dimensional 
methods, and this in spite of the fact that the ionic 
configurations of Na4 and Na^ might let the cylindri- 




Ground state structures for Na|j", 



Na+, 



Na4 



FIG. 2. 

and Nae as obtained in CAPS with the smooth core 
pseudopotential. The numbers report the bond lengths 
in ao, where the uppermost values in each column are 
from the present work, the p*pfj in square brackets from 
Hartree-Fock/CI calculpti|eBS,c3'LZI the ones in parentheses 
from DFT with LDA,Ej'EJ and th c i o ne s in braces from 
all-electron Hartree-Fock calculations.EZlQ 
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cal averaging seem a rather far-fetched approximation. 
But not only the overall geometry, also the bond lengths 
are in good agreement with the ones from the three- 
dimensional calculations. Whereas all-electron Hartree- 
Fock overestimates bond lengths due to missing correla- 
tion effects, the three-dimensional LDA calculations lead 
to the well-known underestimation. By construction, the 
phenomenological pseudopotential compensates this un- 
derestimation, and it is thus reassuring to see that the 
resulting bond lengths are close to the ones found in the 
quantum chemistry calculations. Three-dimensional cal- 
culations—have been performed for a few other neutral 
clusters,B and CAPS is in agreement with the three- 
dimensional geometries in these cases, too. A detailed 
analysis of structures of neutral clusters and their static 
electric polarizability can be found in Ref. |^. 

A further test for structure calculations is obtained by 
comparing the photoabsorption spectra corresponding to 
the theoretically determined structures to the measured 
ones. To this end, we have calculated the valence elec- 
tron excitations for our cluster structures in a collective 
approach. The-ibasic idea of this method has been pre- 
sented earlier ,e3e^ and a detailed discussion of its exten- 
sions and the relation to density funqtipnal theory is the 
subject of a forthcoming publication.Eil Here, we demon- 
strate the validity of our collective model for the test 
cases Na2 and Nas, where experimental data and recent 
ab initio TDLDA calculations are available for compar- 
ison. The first row of plots in Fig. || shows the exci- 
tation spectra obtained with the axial collective model 
for the CAPS structures shown in the insets, the sec- 



Na2 

" Ij 




Nag 






II 

1 _A 






Experiment 


A 


„ J 









3 4 

eV 



FIG. 3. Photoabsorption cross section a in arbitrary units 
against energy in eV for Na2 and Nag. From top to bottom: 
present results for the CAPS structure shown in the inset, 
TDLDA cakulations for the correspopding three-dimensional 
structures,Lj experimental spectra.Lj'Ej The TDLDA and col- 
lective model results have been broadened by 0.06 eV to sim- 
ulate the finite line width of the experiment. 



ond row shows the TDLDA spectra for the correspond- 
ing three-dimensional . gfiop etries, and the bottom row 
shows the experiment .PrP^ The lower two rows were di- 
rectly adapted from Rcf. 
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and the collective model re- 
sults were then plotted on the same scale. Besides small 
differences in the two small peaks seen at high energies 
for Na2, the agreement between the three sets of data is 
very good. This demonstrates that our collective model 
is capable of correctly describing excitations even in small 
systems. It further is to be noticed that the excitation 
energies obtained with the smooth-core pseudopotential 
are « 0.15 eV lower than the TDLDA energies and closer 
to the experimental ones, which can be attributed to the 
larger bond lengths that result with the present pseu- 
dopotential. 



IV. STRUCTURES AND PHOTOABSORPTION 
SPECTRA OF SINGLY-CHARGED CATIONIC 
SODIUM CLUSTERS 

Since the comparisons in the previous section have 
shown that the phenomenological pseudopotential, the 
collective model and CAPS are well suited for the de- 
scription of sodium clusters, a systematic study of clus- 
ters with even electron numbers between 8 and 58 is pre- 
sented. For each cluster a great number of simulated an- 
nealing runs was started from different random configura- 
tions. The search was continued until new runs no longer 
returned new low-energy structures. Although it must 
always be kept in mind that no practical algorithm guar- 
antees that all low-energy structures are found, this ex- 
tensive and unbiased procedure at least gives good hope 
to do so. The aim of this survey is to investigate how elec- 
tronic and ionic shell effects play together in determining 
the cluster structure, and thus gain a better understand- 
ing of how matter grows. The geometry optimization was 
done with CAPS, but the energies of the resulting struc- 
tures have also been obtained in three-dimensional Kohn- 
Sham calculations to check the ordering of isomers with- 
out axial restriction. Photoabsorption spectra are cal- 
culated with the collective model, and comparison with 
the experimental spectra gives further information on the 
relevant structures. 

The "magic" cluster Na.^ is spherical according to 
the jellium model, and the ground state found with 
CAPS is the Ciy structure (a) in Fig. ^. Separated 
by 0.10 eV, CAPS finds the isomer (b), and both 
these geometries were also-Jbund in three-dimensional 
CI and LDA calculations.E3y In addition, CAPS finds 
the third isomer (c) which is higher by 0.23 eV. All 
these clusters have nearly spherical valence electron den- 
sities, and therefore support the shell model picture. But 
the ionic configurations are, of course, non-spherical, 
and this is reflected in the photoabsorption cross sec- 
tions. The first thing to be noted is the overall po- 
sition of the resonance— arb pU | t . wh idi there has beenr-a 
long-standing debate. BEJIlaESEaCaE^ In previous workeS 
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L 







1.5 3.0 3.5 eV 





1.0 1.5 3.0 3.5 eV 1.0 1.5 3.0 

FIG. 4. Low-energy configurations for Na[|", and corre- 
sponding photoabsorption spectra from the collective model. 
A phenomenological Lorentzian line broadening of width 0.08 
eV has been applie^p-iDotted curve: experimental cross sec- 
tion for T=105 KUMi 



we have shown that when ionic structure is taken into 
account via a pseudopotential, the detailed form of the 
potential greatly influences the resonance position. But 
whereas our previous calculations were mainly compared 
to high-temperature data, and, due to limited collec- 
tive model basis sets, fixed the resonance positions only 
within a few percent (as pointed out in Ref . |6|) , we have 
now fully converged the basis. The resulting overall res- 
onance positions are very close to the experimental low- 
temperature positions, showing that consistent inclusion 
of ionic structure gets the plasmon position in the cor- 
rect energy interval. The optical response of . N . alt . has 
also been theoretically investigated previously,ll3c3E3'C3 
and each of these calculations explains the experimental 
spectrum on the basis of a different isomer. In our calcu- 
lations, the spectra of isomers (b) and (c) both reproduce 
the splitting of the main peak seen in the experiment. 
Isomer (a) does not show this feature, but it leads to the 
small sub-peak around 3.5 eV that is also seen experi- 
mentally. The situation becomes transparent when the 
results from our three-dimensional calculations, see Ta- 
ble [II, are taken into account. 



They show that isomer 
(c) is degenerate with isomer (a), and the experimental 
spectrum can thus be explained even at low tempera- 
tures as resulting from a mixture of isomers (a) and (c). 
(The importance of isomerism on dipole spectra in small 
sodium clusters has also been pointed out in Ref. [t^ .) 

The lowest energy structure that we find for Na^^j^ {Did, 
labeled (a) in Fig. |^) can be understood as resulting from 
placing one atom in the center of Naio, or capping iso- 
mer (b) of Na;^ on the quadratic faces. Separated by 
0.11 eV and 0.19 eV, respectively, CAPS finds the 
isomer (b), and its deformed counterpart (c), which is 
only a shallow local minimum in CAPS and easily trans- 
forms into the more stable structure (b). In the three- 
dimensional calculations, the ordering of isomers is re- 
versed: (b) becomes the ground state, and our results 




FIG. 5. Low-energy configurations for Na]*"]^. See text for a 
discussion of photoabsorption spectra. 



are thus consistent with the findings of Ref. ^ We have 
calculated the photoabsorption spectra with the collec- 
tive model and obtain two resonances with heights ap- 
proximately 1:2 for all three geometries. These are cen- 
tered at 2.41/2.94 eV for structure (a), at 2.39/2.94 eV 
for structure (b), and at 2.27/2.98 eV for stmcture (c). 
The experimental photoabsorption spectrumlI3 of Na^j^ 
at 380 K shows two broad peaks with heights 1:2 around 
2.2 eV and 2.85 eV, and a pattern of six peaks when mea- 
sured at 35 K: three low intensity ones at 1.9 eV, 2.2 eV 
and 2.4 eV, and three high intensity ones at 2.5 eV, 2.8 
eV and 3.0 eV. Obviously, the calculated resonances are 
blue shifted by a few percent with respect to the hot ex- 
periment. This is understandable since the calculations 
have been done for T=0 K, and thermal expansion of 
the cluster shifts the plasmon to slightly lower energies. 
(This effect has recently been|piit into evidence quantita- 
tively for the static response.E3) When compared to the 
cold data, the calculated resonances are in the correct 
energy range, but the fine structure that ±lie experiment 
shows is not reproduced. A CI calculatioi£j based on iso- 
mer (b) reproduced some of the experimentally observed 
patterns, but also could not explain all the peajks, and re- 
cent three-dimensional TDLDA calculationalj have lead 
to similar results as the collective model. An explana- 
tion of all details in the experimental spectrum therefore 
might require to consider a mixture of low-energy struc- 
tures, and iMs will be discussed in detail in a separate 
publication.Ea 

For Na^5 CAPS finds four low -energy structures. The 
two lowest ones, (a) and (b) in Fig. |^, are very close in 
energy in both CAPS and the three-dimensional calcula- 
tions, (a) corresponds to a distorted Na^j^ (a) with four 
atoms added along the z-axis, and (b) can be understood 
as a pentagonal bipyramid sharing one ion with Na,^ (a). 
Structure (c) results from structure (a) by moving one of 
the inner single atoms to the bottom face, . and . the oblate 
(d) corresponds to a configuration foundtilO for Nai4, 
but with one atom added along the symmetry axis. In 
agreement with the experimental photoabsorption spec- 
trum, our results point to prolate clusters as the relevant 
structures and thus again verify the prediction of the fde- 
formed, structure averaged jelhum model (SAJM)II3'El 
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3.5 eV 



FIG. 6. Low-energy configurations and photoabsorption 
spectra for Nafg. A phenomenological Lorentzian fine broad- 
ening of width 0.2 eVi4*as been applied. Dotted curve: exper- 
imental cross section.O 



The ground state structure found for Na^^, (a) in Fig. 
^, is close to the one fotuid in the extended Hiickel model 
for the neutral cluster It consists of an icosahedral core 
with a crown of four atoms, and the only difference to the 
result of Ref. ^ is that the cylindrical averaging forces 
the crown atoms into a square. This shows that for such 
sizes the differences between charged and neutral clus- 
ters can already be small. Structure (b), which is slightly 
higher in energy than (a) in the three-dimensional calcu- 
lations, has an even more prolate valence electron den- 
sity, and the collective resonances for these clusters lie at 
2.59 eV and 2.86 eV [structure (a)] and at 2.34 eV and 
2.95 eV [structure (b)]. They are thus in contrast to the 
experimental spectrum that points to an oblate cluster. 
CAPS finds one more isomer, structure (c), that indeed 
leads to an oblate valence electron density. In CAPS, it 




1.75 2.50 115 3.00 3.75 GV 



FIG. 7. Low-energy configurations for Najt,, and pho- 
toabsorption spectrum for isomer (c). A phenomenological 
Lorentzian line broadening of width 0.22 ejV, has been applied. 
Dotted curve: experimental cross section.lIj 



is 0.16 eV higher than (a), with a half occupied orbital. 
The collective spectrum for this isomer is very close to 
the experimental one, as shown in Fig. 0. However, also 
the three-dimensional calculation gives a non-negligible 
energ y di fference between (c) and the prolate isomers, cf. 



Table III. Therefore, in this case a three-dimensional re- 
laxation of the ions would be necessary to check whether 
Jahn- Teller distortions lower the energy of isomer (c) so 
much that it can account for the experimental spectrum. 

For Naj'g CAPS finds the four stable geometries shown 
in Fig. ||. The double-icosahedron (a) is the ground state, 
and it is separated from structure (b) by 0.23 eV, from 
(c) by 0.30 eV, and from (d) by 0.35 eV. However, in 
the calculations without axial restriction, the energetic 
differences are considerably reduced: isomer (d) becomes 
equivalent to isomer (b), and both are separated from the 
ground state by only 0.12 eV. The experimental photoab- 
sorption spectrum sheds further light on the situation. It 
shows a high peak at about 2.7 eV, followed by a lower 
one at about 2.9 eV, and is thus very similar to the col- 
lective spectrum of isomer (d). The spectra of the other 
isomers all give the peaks in reversed order and thus do 
not resemble the experiment. (The collective spectrum 
of (c), which is not shown in order not to overload Fig. 

shows a small peak at 2.84 eV and a higher one at 
2.93 eV.) Since an internal cluster temperature of about 
60 K is sufficient to go from isomer (a) to (d), and since 
the temperature in the experiment is about 105 K, it is 




FIG. 8. Low-energy configurations and photoabsorption 
spectra for Najj. A phenomenological Lorentzian line broad- 
ening of width 0.15 eV ha^ been applied. Dotted curve: ex- 
perimental cross section.Lj 
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plausible that isomer (d) will contribute to the spectrum. 
The experimental spectrum however looks as if it is domi- 
nated by isomer (d). One explanation for this can be that 
structure (d) has a larger "catchment area" than (a) and 
is entropically favored. The second possibility, which is 
suggested by the relatively large difference between the 
energy found in CAPS and the three-dimensional calcu- 
lation, is that a three-dimensional relaxation of the ions 
would turn structure (d) into the ground state. 



3.5 eV 




3JeV 



3.5 eV 



FIG. 9. Low-energy configurations and photoabsorption 
spectra for Najj^. A phenomenological Lorentzian fine broad- 
ening of width 0.15 eV Ips been applied. Dotted curve: ex- 
perimental cross section.Lj 

Two degenerate structures are found for NaJ,^ in 
CAPS, and also the three-dimensional calculation leads 
to nearly equal energies. Structure (a) results from the 
double icosahedron of Na^g (a) by adding one ion to each 
of the lower tow pentagons, and (b) results from (a) by 
moving the lowest ion to the uppermost pentagon. Both 
geometries lead to collective spectra with two main tran- 
sitions and some strength above the main reaauances. 
This is in agreement with the experimental dataOO that 
also show this high-energy tail. Concerning the main 
transitions, the two isomers are somewhat different: in 
structure (a) the main peaks are separated by about 0.1 
eV, whereas they nearly fall together for structure (b). 
Since the structures are extremely close in total energy 
and since the main peaks of (b) energetically fall exactly 
between the main peaks of (a), the experimental spec- 
trum can consistently be explained as resulting from a 
mixture of both isomers. 

The simulations performed for Najg led to strongly 
prolate shapes. Th^ is partially understandable on the 
basis of the SAJMiLl the valence electrons force the clus- 
ter into the prolate regime. But the inclusion of ionic 
structure even enhances the prolate tendency, because 
the lowest configuration, shown in Fig. |l^, is a triple 
icosahedron, and in order to build this configuration, a 
quadrupole moment larger than the jellium prediction is 
necessary. From this observation it also becomes clear 




2,25 2.50 2.75 3.00 3,25 3.50 



FIG. 10. CAPS ground-state configuration and photoab- 
sorption spectrum for Na^^g. A phenomenological Lorentzian 
line broadening of width 0.25 eVpbas been applied. Dotted 
curve: experimental cross section.Lj 



that the pentagonal bipyramid which has bees jobserved 
as a building block of the smaller clustcrsllaO'Ell is also 
important for the larger sizes. For the triple icosahe- 
dron, collective resonances are obsepjed at 2.44 eV and 
2.96 eV, and also the experimentsOiIj indicate peaks at 
these energies. But the latter shows an additional peak 
around 2.75 eV which is not found in the calculations. 
Comparing this to the situation encountereccS for Na27 
shows that the geometries of the two clusters are closely 
related: Na2'7 triple icosahedron with two atoms 

added to the central pentagons. This shows the con- 
sistency of the structure calculations, and for Na^y the 
collective spectrum nicely matches with the experimen- 
tal data. The triple icosahedron is energetically strongly 
favored over the other structures that CAPS finds. All 
of them are variations of the ground state geometry with 
local distortions, and a three-dimensional calculation was 
performed for the structure that in CAPS is the second 
lowest. Here, the difference in total energy is a little 
smaller than in CAPS, but still large, cf. Table 111. Since 
furthermore all low-energy structures give rise to collec- 
tive spectra that are similar to the one shown in Fig. nO, 
the explanation for the middle peak observed for NaJs 
remains an open question. 

Comparing the structures of Na^^ and Naj y to the 
ones found for Na+ , NaJ^ and Na+ (Fig. [uf Fig. |l|, 
and Refs. ^,^6|) allows to identify a growth pattern: the 
triple icosahedron serves as a building block for the larger 
clusters. Additional atoms are added on outside faces, as 
seen in the ground states of Naa'y and Na43, or can be 
packed between the pentagonal subunits, as observed for 



the isomers of Najy, where the two additional ions are 
placed at different positions "inside" the Na^g structure. 
The collective spectra of the low-energy structures (a) - 
(c) of NaJ]^ reflect the prolate electron densities of these 
clusters and show one peak around 2.6 eV, carrying about 
30% oscillator strength, another one around 2.95 eV car- 
rying about 50% strength, and the rest of the strength 
scattered at higher energies. Structure (d) is not stable 
against Jahn Teller distortions in three dimensions. In 



9 




a b c d 

FIG. 11. Low-energy configurations for NaJ^. 



view of the cluster size, the CAPS differences in total en- 
ergy for Na^]^ are rather small, and they are even smaller 
in the three-dimensional approach. Since the temper- 
ature in the measurement of the photoabsorption cross 
section was again at least 105 K, it can be concluded from 
the calculations that a variety of isomers will contribute 
to the experimental spectrum. It is t]|ijis,not astonishing 
that the available experimental dataOO do not resolve 
separated peaks but show a rather broad bump. 




2.25 1.50 1J5 3.00 3.15 eV 



FIG. 12. CAPS ground state configuration and photoab- 
sorption spectrum for Najg. A phenomenological Lorentzian 
line broadening of width 0.15 eVrias been applied. Dotted 
curve: experimental cross section." 

The pronounced deformation seen for Na43 ^i§- |H 
results from an interplay between ions and valence elec- 
trons. The electrons "push" the cluster into the octupole, 
i.e. pear-shaped form, and the ions arrange under this 
"constraint" . But the ions favor the icosahedral core, and 
this increases the deformation. The octupole moment 
therefore is larger by a factor 1.15 than in the SAJM.LlI 
The photoabsorption spectrum calculated for Na^g shows 
a strong peak at 2.7 eV, followed by a high-energy tail. It 
is in close agreement with the experimental data, which 
thus support our structure calculations. 

Based on the shell and deformed jellium model a con- 
figuration would be expected for Na^g that gives rise to a 
prolate valence electron density. The CAPS calculations, 
however, consistently led to nearly spherical or slightly 
oblate clusters as low-energy configurations. The low- 
est energy was found for structure (a) in Fig. 0[ It 




L13 1.50 1.75 3.00 3.15 3.50eV 



FIG. 13. Low-energy configurations and corresponding 
photoabsorption spectra for Na^g. A phenomenological 
Lorentzian line broadening of width 0.17 eV hcpibeen ap- 
plied. Dotted curve: experimental cross section.LJ See text 



for a discussion of further experiments! 



has a fivefold symmetry axis and again shows the close- 
to-icosahedral core discussed previously. Within 0.4 eV 
structures (b) and (c) are found with valence electron 
densities that are more oblate. These structures show 
half-occupied highest orbitals, indicating that they would 
Jahn- Teller relax if the axial restriction on the electrons 
were lifted. The lowest prolate isomer that was found, 
structure (d), has a quadrupolejjuoment that is close to 
the one predicted by the SAJMQ but in CAPS this con- 
figuration is 0.73 eV higher than the ground state. The 
annealing was also started from an icosahedron with the 
nearest-neighbor distance of bulk sodium. In this cal- 
culation, the strictly geometrical icosahedron transforms 
into structure (a), which can be regarded as a slightly 
distorted icosahedron. The experimental photoabsorp- 
tion spectrum of Na^g as measured by the Freiburg group 
is shown in the lowest panel of Fig. and the. spec- 
trum has also been measured by Meibom et al£3 Both 
spectra have in common that one broad peak is seen, fol- 
lowed by a second, smaller peak or a high-energy shoul- 
der. This is in contrast to the prediction of the jellium 
model, because the prolate structures found there give 
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the peaks in reversed order, similar to the spectrum of 
isomer (d). The CAPS low-energy structures, however, 
lead to collective spectra with a high peak followed by 
a lower one, and are thus in better agreement with the 
experiment. Tiie total oscillator strength measured in 
the experimentEa was 70% - 80%. This also agrees with 
the collective model calculations that find 79% strength 
within the plotted range for structure (a), 76% for (b) 
and 80% for (c). As a test for the collective model, a 
TDLDA calculation with excitation in z-direction was 
performed for structure (a). In the TDLDA, the high- 
energy shoulder is more pronounced than in the collective 
model. To further clear up the situation, the structures 
(a), (b) and (d) were relcised in Born-Oppenheimer ab 
initio molecular dynamics JIj i.e., fully threej-dimensional 
with the TrouUier-Martins pseudopotentialE^ The bond 
lengths of all three structures shrinked by a few percent 
due to the different pseudopotential, and isomers (b) and 
(d) distorted slightly, but otherwise the geometries stayed 
unaltered. As seen in Table III, the differences in total 



energy are somewhat smaller in the ab initio calculations, 
but the ordering of isomers is the same as in CAPS. Fi- 
nally, as reported in Ref. |2^, the CAPS structure (a) was 
annealed for 10 ps at about 220 K. In this annealing, the 
cluster became even more similar to an icosahedron and 
its overall shape oscillated around the nearly spherical 
shape of the icosahedron. Thus, the ab initio calcula- 
tions confirm the finding that in contrast to the deformed 
jellium model prediction, Na^Jg at low to moderate tem- 
peratures has a close to spherical valence electron density. 

The CAPS results for Na^y, see Fig. |l^, are consistent 
with the results for Na^5 : again the low-energy structures 
are not prolate, but spherical or slightly oblate. The col- 
lective spectra somewhat underestimate the high-energy 
shoulder, and this can be attributed to selective particle- 
hole effects just as in the case-|p£ Na^5. But the over- 
all agreement with experimentllallj is considerably better 
than in the deformed-jellium calculation. Structure (a) 
has the lowest energy, followed by isomers (b) and (c) 
that are higher by 0.22 eV. For structure (c) we observe 
a half-occupied orbital. With respect to the ionic geom- 
etry, structure (a) is similar to the third isomer of Na^g , 
and (c) corresponds to the ground state of Na^5 with the 
two additional ions added on top and bottom faces. An- 
other isomer, not shown in Fig. |lj for the sake of clarity, 
is found, and there the two additional ions are added in 
the equatorial plane of Nagj (a). 

In the case of Na^g the search for low-energy structures 
has not been as exhaustive as for the other clusters. But 
several annealing runs were started from random config- 
urations, and the low-energy geometries that were found 
again showed the ions arranged rather regularly and sim- 
ilar to the structures just discussed. Therefore, further 
simulations were started from the geometries found for 
Najjg and Na^y, plus four or two atoms, respectively, that 
were added at randomly chosen sites. One of these runs 
led to a geometry that corresponds to Na^^ (a), capped 
by two atoms on top and bottom. This structure has the 
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FIG. 14. Low-energy configurations and corresponding 
photoabsorption spectra for Nag^. A phenomenological 
Lorentzian line broadening of width 0.17 eV hcp-ibeen ap- 
plied. Dotted curve: experimental cross section.Cj See text 
for a discussion of further experimentsES 



lowest energy of all that were found, and it does not show 
signs of Jahn-TcUer instability. The comparison between 
the experimentalLj and theoretical photoabsorption spec- 
tra fits into the previous discussion. The strongest collec- 
tive resonance is seen at 2.8 eV, followed by a smaller one 
at 2.9 eV and a little strength scattered around 3.25 eV, 
i.e., the collective model leads to qualitative agreement 
with the experimental data. 



V. DISCUSSION AND CONCLUSION 

The systematic survey showed that for most of the 
smaller Na clusters, the overall deformation is determined 
by electronic shell effects even when the ionic structure 
is explicitly included. This explains the great success of 
the deformed jellium model. However, the survey at the 
same time clearly exhibits the limitations of the jellium 
picture. Besides the fact that some spectroscopic pat- 
terns, e.g., as seen for Na^, are directly related to details 
in the ionic configuration, the fundamental improvement 
brought about by taking into account the ionic struc- 
ture is that growth systematics can be identified. In 
the smallest clusters, the pentagonal bipyramid is a fre- 
quently appearing structure and is seen, e.g., in Na^", 
Na^g, NaJ^. In Nafg, three of these smallest units are 
combined to build up the double icosahedron, and from 
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then on, the electronic shell effects and the preference of 
the ionic structure for icosahedral geometries work hand 
in hand in determining the cluster structure. This is seen 
most impressively for Na43. ^^55^ however, the situa- 
tion changes. Due to the ionic shell closing, the influence 
of the ions wins over the electronic shell eflfects, resulting 
in nearly spherical structures. This finding is supported 
by the experimental photoabsorption data and has been 
verifled in ab initio calculations. 

In summary, we have presented a local pseudopoten- 
tial for sodium that accurately models the core-valence 
interaction and correctly describes the atom, the bulk, 
and finite clusters. We have demonstrated the influence 
of pseudopotentials on structural properties and the di- 
rect influence of the bond lengths on the resonance po- 
sitions. This shows that even for the most simple metal, 
a pseudopotential must be constructed carefully. Cluster 
structures were calculated in axially averaged and three- 
dimensional Kohn-Sham calculations. A detailed com- 
parison with ah initio work demonstrated the validity 
of the CAPS as a tool for the approximate determina- 
tion of cluster configurations. The systematical survey 
for clusters with up to 59 ions revealed an icosahedral 
growth pattern. Collective resonances were calculated 
for the theoretically determined structures, and compar- 
ison with the experimental photoabsorption spectra con- 
firmed the results of our structure calculations. Thus, a 
distinct step in the growing process of matter has been 
put into evidence, namely the transition from electroni- 
cally to ionically determined configurations. 
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Cluster 


Ai?pAPs/eV 


AE-m/eV 








Na;^ b 


0.10 


0.12 


Nag c 


0.23 


G' 


IN a^ 2. a 




U.lu 


Na+ h 


U. J- J- 


n 

Vjr 


NaJ^ c 


0.19 


0.10 


Na+ a 


G 


0.04 


NaJ; b 


0.03 


G 


Na+ c 


0.15 


0.16 


Na+ d 


0.23 


0.15 


Na^7 a 


G 


G 


Na+ b 


G' 


0.08 


NaJ^ c 


0.16 


0.20 


Na^ a 


G 


G 


Na+ b 


0.23 


0.12 


Na+g c 


0.30 


0.20 


Na+g d 


0.35 


0.12 


Najj^ a 


G 


G 


Na+ b 


G' 


0.03 


Najg a 


G 


G 


Na+ b 


0.41 


0.29 


Najj a 


G 


G 


NaJi b 


0.14 


0.08 


NaJi c 


0.19 


0.10 


Na+ a 


G 


G 


Na+ b 


0.41 


0.30 


Na+ d 


0.73 


0.31 



TABLE III. Differences in total energy for NaJ clusters. 
Small letters behind the cluster symbol label structures and 
refer to Fig. ^ - Fig. |l^ and the main text. G (and G') de- 
note the structure with lowest energy for a given cluster size. 
The left column for each size gives the difference in total en- 
ergy between the corresponding structure and G as found in 
CAPS. The right column gives the energetic differences found 
for the same structures in the three-dimensional Kohn-Sham 
calculation. The 3D values for Na^g were obtained by re- 
laxing the CAPS stpHctures in ab initio Born-Oppenheimer 
molecular dynamics,EJ see text for discussion. 
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